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Numerical approximation of the 
Euler-Poisson-Boltzmann model 
in the quasineutral limit 
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Abstract 

This paper analyzes various schemes for the Euler-Poisson-Boltzmann (EPB) 
model of plasma physics. This model consists of the pressureless gas dynamics 
equations coupled with the Poisson equation and where the Boltzmann relation 
relates the potential to the electron density. If the quasi-neutral assumption is 
made, the Poisson equation is replaced by the constraint of zero local charge and 
the model reduces to the Isothermal Compressible Euler (ICE) model. We compare 
a numerical strategy based on the EPB model to a strategy using a reformulation 
(called REPB formulation). The REPB scheme captures the quasi-neutral limit 
more accurately. 
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1 Introduction 



The goal of this paper is to analyze various schemes for the Euler-Poisson-Boltzmann 
model of plasma physics. The Euler-Poisson-Boltzmann (EPB) model describes the 
plasma ions through a system of pressureless gas dynamics equations subjected to an 
electrostatic force. The electrostatic potential is related to the ion and electron densities 
through the Poisson equation. The Boltzmann relation provides a non-linear relationship 
between the potential and the electron density which allows to close the system. More 
precisely, the Euler-Poisson-Boltzmann model is written 



Here, n(x,t) > 0, u(x,t) £ M. d and 4>{x,t) e R stand for the ion density, ion velocity and 
electric potential respectively, which depend on the space- variable x G M. d and on the time 
t > 0. We suppose that the ions bear a single positive elementary charge e and we denote 
by m, their mass. The electron temperature T is supposed uniform and constant in time, 
eo and ks respectively refer to the vacuum permittivity and the Boltzmann constant. The 
operators V, V- and A are respectively the gradient, divergence and Laplace operators 
and u £g> u denotes the tensor product of the vector u with itself, n* is usually fixed by 
either imposing zero total charge 



or by assuming that the net charge is zero at a given point x* (for instance at the bound- 
ary): 



For simplicity, this work is restricted to dimension d = 1 but the concepts extend to di- 
mensions d > 2 without additional difficulties and numerical applications will be reported 
in future work. 

The ion pressure force is neglected. This is a commonly made assumption in plasma 
physics [5j [20] for elementary text books of plasma physics. The inclusion of an ion 
pressure term would not modify the subsequent analysis and is omitted for simplicity. 
Additionally, the pressureless system has interesting multi-valued solutions which are lost 
in the case of the non-pressureless model [231 12^1 125] . 

If the quasi-neutral assumption is made, the Poisson equation (11.31) is replaced by the 
constraint of zero local charge: 





(1.2) 






In this context, we can write 



* ( ^ \V7A kBT Y7( * I ^ \\ kBT V7 ft A\ 

nV0 = n exp — — V0 = V n exp — — = Vn, 1.4 

k B l e k B l e 

and the quasi-neutral Euler-Poisson-Boltzmann model coincides with Isothermal Com- 
pressible Euler (ICE) model: 

d t n + V • (nu) = 0, 

m(d t (nu) + V • (nu ® u)) + V(nk B T) = 0. 

The passage from EPB to ICE can be understood by a suitable scaling of the model, 
which highlights the role of the scaled Debye length: 



Ad > / e k B T 

~T: A D — I — 5— 



where L is the typical size of the system under consideration. measures the spatial 
scale associated with the electrostatic interaction between the particles. The dimensionless 
parameter A is usually small, which formalizes the fact that the electrostatic interaction 
occurs at spatial scales which are much smaller than the usual scales of interest. However, 
there are situations, for instance in boundary layers, or at the plasma-vacuum interface, 
where the electrostatic interaction scale must be taken into account. This means that the 
choice of the relevant macroscopic length L may depend on the location inside the system 
and that in general, the parameter A may vary by orders of magnitude from one part of 
the domain to another one. The scaling will be presented in more detail in section [2j 

This paper is concerned with discretization methods for the EPB model in situations 
where A can vary from order one to very small values. Therefore, the targeted schemes 
must correctly capture the transition from the EPB to the ICE models. With this ob- 
jective in mind, we will compare two strategies: a first one which uses the EPB model 
in its original form, and a second one which reformulates the EPB model in such a way 
that it explicitly appears as a perturbation of the ICE model. This reformulation uses 
the Poisson equation in the form: 

n = n exp(— — ) A0. 

k B T e 

Then, with the same algebra as for (11.41) . we get: 

nV(p = n*exp(-^)V0- — 
k B T e 

= ^^))-V-w«w)-v(K)) 

e k B l e 2 

= ^Vn + ^VA0-^(V-(V0®V0)-V(^)). (1.6) 

e e z e 2 
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We note that this expression is reminiscent of the Maxwell stress tensor. Then, the 
EPB model is equivalently written as the following reformulated Euler-Poisson-Boltzmann 
(REPB) model: 

d t n + V • (nu) = 0, (1.7) 
m(d t (nu) + V • (nu <g> u)) + V(nk B T) = 

= - eo V(^A0 + + e V • (V0 ® V0), (1.8) 

-A0 = -(n-n*exp(-^)). (1.9) 
e k B T 

In this way, the ICE appears at the left-hand side of ( 11.7j) . (11.81) . The scaling analysis 
will show that the right-hand side of (11.81) is of order A 2 and is therefore negligible (if the 
gradients of the potential are smooth) in the limit A — > 0. 

The goal of this paper is to propose and analyze two schemes for the EPB model which 
provide the correct ICE limit when A — > 0: the first one is based on the initial formulation 
EPB and the second one, on the reformulated form REPB. To present the schemes, we 
note that both EPB and REPB can be put under the form 

d t n + V ■ (nu) = 0, 

m(d t (nu) + V ■ (nu <g> u)) + Vp(n) = S(n, V0), 

-A0= — (n-n*exp(-^— )), 
e k B l 

with 

p(n) = 0, S(n, V0) = — enV(p, 

in the case of EPB and 

p(n) = nk B T, S(n, V0) = -e o V(^A0 + i^E) + e V • (V0 <g> V0), 

in the case of the REPB. 

Both schemes use the following time-semi-discretization which is implicit in the Poisson 
equation and in the source terms of the momentum equation: 

5~\n k+1 -n k ) + V- ((nu) k ) = 0, 

m(8- l ((nu) k+1 - (nuf) + V ■ ((nu) k ® u k )) + Vp(n k ) = S(n k+1 , V/ +1 ), 

-A/ +1 = -(n fc+1 -n*exp(^)). 
e k b 1 

Here, 5 is the time step and the exponent k e N refers to the approximation at time 
t k = k5 (i.e. n k (x) ~ n(x,t k ), ...). For the EPB model, this discretization is classical 
[T8] . For both systems, in spite of its implicit character, the recursion can be solved in 
an explicit way, by first updating the mass equation to find n k+1 , then using the Poisson 
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equation to find the potential <p k+1 and finally using the momentum equation to find u k+l . 
Of course, the space operators have to be discretized as well and a simple shock-capturing 
method is used, namely the Local Lax-Friedrichs or Rusanov scheme [221 ITS] . A time- 
splitting is used where the hydrodynamic equations are evolved without the source terms. 
Then, the Poisson equation is solved with the value of the density found at the end of the 
first step of the splitting. With the newly computed value of the potential, the evolution 
of the hydrodynamic quantities due to the source terms are computed. 

In section[3l we will show that both schemes are Asymptotic-Preserving. The Asympto- 
tic-Preserving property can be defined as follows. Consider a singular perturbation prob- 
lem P x whose solutions converge to those of a limit problem P° when A — > (here P x is 
the EPB model and P° is the ICE model). A scheme P$ h for problem P x with time-step 
S and space-step h is called Asymptotic Preserving (or AP) if it is stable independently of 
the value of A when A — > and if the scheme P$ h obtained by letting A — >■ in P$ h with 
fixed (5, h) is consistent with problem P°. This property is illustrated by the commutative 
diagram below: 

pA (*.fr)-+0 pA 

A^O A^O 
pO (W^O p0 

The concept of an AP scheme has been introduced by S. Jin [TH] for diffusive limits of 
kinetic models and has been widely expanded since then. 

In section [3j we will perform a linearized stability analysis which shows that both 
schemes are stable independently of A in the limit 5 — > 0, provided that the usual CFL 
condition of the ICE model is satisfied. However, the scheme based on the REPB for- 
mulation has several advantages over the one based on the EPB form. A first advantage 
lies in the fact that the hydrodynamic part of the EPB model is a pressureless gas dy- 
namics model, which is weakly unstable and can produce delta concentrations (see e.g. 
El IH El [7]). By contrast, the hydrodynamic part of the REPB is the usual ICE model, 
which is strongly hyperbolic, and which is much more stable than the pressureless gas 
dynamics model. A second advantage is the fact that the limit A — > of the REPB- 
based scheme provides a conservative discretization of the ICE model, while that of the 
EPB-based scheme leads to a scheme in non- conservative form. We can expect that the 
accuracy of the PB-based scheme degrades when A 1 for solutions involving disconti- 
nuities, which is not the case for the REPB-based scheme. 

To illustrate the theoretical findings, one-dimensional numerical simulations are pre- 
sented in section following a discussion of the spatial discretization in section |H First, 
setting A = 1, analytic solutions can be derived in the form of solitary waves thanks to the 
Sagdeev potential theory [5]. Both schemes are compared to these analytical solutions, 
and show a similar behavior, with a slightly larger numerical diffusion in the case of the 
REPB scheme. Then a Riemann problem test case consisting of two outgoing shock waves 
is investigated. In the A 1 regime, the REPB scheme captures the right hydrodynamic 
shocks while the EPB scheme develops spurious oscillations. This confirms the better 
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behavior of the REPB scheme in the A < 1 regime. Finally, a test-problem related to 
multivalued solutions and proposed in [231 [21] is investigated. In this case, both schemes 
shows a similar behavior. To summarize, the REPB scheme captures well the A 1 limit 
but is slightly more diffusive in the A = 0(1) regime. But the extra numerical diffusion is 
mild and then shows that the REPB scheme is superior to and should be preferred over 
the EPB scheme in most situations. 

We conclude this section by some bibliographical remarks. The Euler-Poisson-Boltzmann 
model has been recently analyzed in the context of sheath dynamics [22] arid multi- valued 
solutions have been computed using level set methods [231 121] • Numerical schemes for the 
quasineutral limit of plasma problems has been the subject of vast literature, mostly for 
the Vlasov-Poisson equation and for Particle-in-Cell (PIC) methods. It is virtually im- 
possible to cite all relevant references and we only refer to the seminal ones [U [211 I2S1 [27] • 
For fluid models of plasmas, the literature is comparatively less abundant. We can refer 
to the pioneering work [12] , and more recently to [SI QUI [23 ED] • Recently, AP-schemes for 
the two- fluid Euler-Poisson model [TT] [T5] or Vlasov-Poisson model [U [T2] in the quasineu- 
tral limit have been proposed. The drift-fluid limit of magnetized plasmas has also been 
considered [131 EH] as well as other applications such as small Mach-number flows [TT] . 
However, none of these works is concerned with Boltzmannian electrons. 



In this section we present a derivation of the EPB model from the two fluid Euler-Poisson 
system and we introduce its scaling. From now on, we will restrict ourselves to one- 
dimensional models. 

2.1 Derivation of the EPB model 

We consider a plasma composed of two species of charged particles: positively charged 
ions and electrons. The ions are supposed singly charged. The modeling of such a plasma 
by means of fluid equations uses a system of compressible Euler equations for each species, 
coupled by the Poisson equation. We denote by m i>e the ion and electron masses, n; >e the 
ion and electron densities and Ui >e the ion and electrons mean velocities. We denote by e 
the elementary charge (i.e. the ion charge is +e > and the electron charge is — e < 0). 
We assume that the ion temperature is negligible so that the ions follow a pressureless gas 
dynamics model. Electrons are assumed isothermal with a non-zero constant and uniform 
temperature T e . Then the electron pressure law satisfies p e = n e kBT e where k% denotes 
the Boltzmann constant. The balance laws for both species are given by 



2 The EPB model and the scaling 



d t rii + d x (riiUi) = 0, 

mi(d t {niUi) + d x {riiul)) = -enid x <f>, 

d t n e + d x (n e u e ) = 0, 

m e (d t {n e u e ) + d x {n e u 2 e )) + d x {n e k B T e ) = en e d x 4>, 



(2.1) 
(2.2) 
(2.3) 
(2.4) 
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where <fi denotes the electric potential. It satisfies the Poisson equation: 

- e d%<p = e(rij - n e ), (2.5) 

where e is the vacuum permittivity. 

The electrons being much lighter than the ions, it is legitimate to take the limit m e — > 
in the electron momentum equation. In this limit, we formally obtain: 

d x (n e k B T e ) = en e d x (f>. 

Integration with respect to x leads to the Boltzmann relation: 

n e = n*exp(-^), (2.6) 

where n* is fixed by some condition (e.g. vanishing total charge or vanishing local charge 
at one given point such as a boundary point, see section [T|). The Boltzmann relation 
shows that the electron density automatically adjusts to the potential. The EPB model 
therefore consists of the ion mass and momentum conservation equations (12. ip . ( 12. 21) . 
the Poisson equation ( 12. 5ft and the Boltzmann relation ( 12.61) . With a change of notation 
rii — > n, Ui — > u, mi — > m, T e — > T, we find the EPB model (ll.ll) - (ll.3p which has been 
introduced in section [TJ 

We note that, in the present one-dimensional setting, the electron velocity u e can 
be computed from the electron density equation (12. 3p . thanks to the value of n e and 
therefore, of <f), obtained by the resolution of the EPB model. Therefore, the computation 
of u e is decoupled from the computation of the other unknowns n«, Ui and involved in 
the EPB model and will be discarded in the present work. In two or higher dimensions, 
the computation of u e requires the resolution of the electron momentum equation, which 
takes the form in the small electron mass limit: 

d t {n e u e ) + V ■ (n e u e ® u e ) + n e Wip = 0, (2.7) 

where ip = lim me _ i ,o( m ^ 1 (^B^ hi n e — e<f))). The quantity ip plays the same role as the 
pressure in the incompressible Euler equation. It is computed thanks to the electron 
density equation (12. 3p which appears as a (non-zero) divergence constraint on u e . In 
this sense, the limit m e — > is similar to the 'low Mach-number' limit of isentropic 
compressible gas dynamics. However, the question of the resolution of (12 .7p is left to 
future work. 

2.2 Scaling of the EPB model 

In this section, we return to the EPB model in the form (ll.ip - (ll.3l) and, with the notations 
of section CD, we introduce a scaling of the physical quantities. Let xo, t , u , O and n be 
space, time, velocity, potential and density scales. Scaled position, time, velocity, potential 
and density are defined by x = x/xq, i = t/to, u = u/u , = — 4>/4>o and n = n/n . 
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We choose xq to be the typical size of the system (for instance an inter-electrode distance 
or the size of the vacuum chamber). The velocity scale is chosen equal to the ion sound 
speed Mo = (fceT/m) 1//2 . We note that the ion sound speed is constructed with the ion 
mass but with the electron temperature. This is clear from the ICE model (see section 
[p. We also choose no = n*. Finally, <po = ksT/e is the so-called thermal potential. Note 
that we have introduced a sign change in the potential scaling because we find it more 
convenient to work in terms of the electron potential energy rather than in terms of the 
electric potential. 

Inserting this scaling and omitting the bars gives rise to the EPB model in scaled 
form: 

d t n x + d x (n x u x ) = 0, (2.8) 
d t (n x u x ) + d x {n x u x u x ) = n x d x (f) X , (2.9) 
\ 2 dl<p x = n x -e-' t >\ (2.10) 

where A is the scaled Debye length (jl.5p . 

It will be useful to consider the linearized EPB model about the state defined by 
n x = 1, u x = 0, <p x = (which is obviously a stationary solution). Expanding n x = l+en x , 
u x = eu x and <p x = e<f) x , with e C 1 being the intensity of the perturbation to the 
stationary state, and retaining only the linear terms in e, we find the linearized EPB 
model: 

d t n x + d x u x = 0, (2.11) 
d t u x = dj x , (2.12) 
\ 2 d 2 J x = n x + 4> x . (2.13) 

Introducing n x , u x , (f) x , the partial Fourier transforms of h x , u x , (f> x with respect to x, we 
are led to the system of ODE's: 

d t n x + i£u x = 0, (2.14) 
d t u x = i^ x , (2.15) 
-A 2 «e 2 A = n x + 0\ (2.16) 

where £ is the Fourier dual variable to x. We note that the general solution of this model 
takes the form 

h x \ n } _ 



ii x I \ u'_ 



f • (2.17) 



with 



4 = ±77-^Lw^ (2-18) 



(1 + A 2 £ 2 )V2 ; 

and n^, are given functions of £, fixed by the initial conditions of the problem. In 

.A 



particular, since s+ are pure imaginary numbers, the L 2 norm of the solution is preserved 



with time. 

Now, we investigate the quasi-neutral limit A — > in the next section. 
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2.3 The quasineutral limit: the ICE model 

Formally passing to the limit A — > in the EPB model in scaled form and supposing that 
n \ _^ n o^ u \ _^ u o^ ^o^ we are ^ e( j ^ Q following model: 

d t n° + d x (n°u°) =0, 

d t (n°u°) + d x (n°u°u°) = n°d x 0°, 

As a consequence of the last relation (which imposes to the ions to satisfy the Boltzmann 
relation of the electrons), we can write (see also f ll.4j) ): 

n°d x( j> = e-^0° = -d x (e^°) = -d x n°, (2.19) 

and, inserting this relation into the momentum equation leads to: 

d t (n°u°) + d x (n°u°u°) + d x n° = 0. 

Therefore, the quasineutral limit A — > consists of the Isothermal Compressible Euler 
system (ICE) complemented by the Boltzmann relation for the potential: 

d t n° + d x {n°u°) = 0, 

d t (n°u°) + d x (n°u°u°) + d x n° = 0, 

n° = e-*°. 

Similarly to the EPB model, the ICE model can be linearized about the state defined 
by n° = 1, u° = 0, 0° = 0. We find (with the same notations as for the EPB model), in 
Fourier space: 

d t n° + i£u° = 0, 
d t u° + i^h° = 0. 

The general solution of this model takes the same form (12.171) as for the linearized EPB 
model with s^. replaced by = ±z£. We note that (see (12.181) ) 

s° + = lim si. 

Therefore, the wave speeds of the linearized EPB model converge to those of the linearized 
ICE model (which are nothing but the acoustic wave speeds). There is no singularity of 
the limit A — > as regards the wave-speeds. In this respect the quasineutral limit A — > 
is not a singular limit. This fact contrasts with the situation of the quasineutral limit 
of the 2-fluid Euler system, where the electron plasma oscillation frequency converges to 
infinity. Therefore, we expect that the numerical treatment of the quasineutral limit in 
the Euler-Poisson-Boltzmann case will be easier. 
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2.4 Reformulation of the EPB model 

To better capture the transition from the EPB model to the ICE model, it is useful to 
reformulate the EPB model in such a way that it explicitly appears as a perturbation of 
the ICE model. Using (I2.10p . in the spirit of f !2. lQj) . we can write (see also ( II. 6p ): 

n'd^ = (e^ A + A 2 A ) d x <j> x 

= 9,(-e^ + y(40 A ) 2 ) 

= ^(A 2 A -n A + y(^) 2 )- 

We note that some simplification arises in the 1-dimensional case, compared to the multi- 
dimensional case of ( II. 6p . Inserting this expression in the momentum equation leads to 
the reformulated EPB systems (REPB): 

d t n x + 4(nV) = 0, (2.20) 

d t (n x u x ) + d x (n x u x u x ) + d x n x = X 2 d x (d 2 ^ + ^0 A ) 2 ) , (2-21) 

n x - = \ 2 d 2 x (j) x . (2.22) 

In this formulation, the ICE model explicitly appears at the left-hand side of (I2.20p - 
(12.221) . Additionally, the remaining terms, at the right-hand side of the equations are 
formally of order A 2 . Therefore, the EPB model explicitly appears as an order 0(A 2 ) 
perturbation of the ICE model. 

We stress the fact that the REPB model is equivalent to the original EPB model. 
However, at the discrete level, schemes based on the REPB model may differ from those 
based on the EPB model. The goal of the present article is to compare the properties 
of schemes based on these two formulations, in relation to their Asymptotic-Preserving 
(AP) properties when A — > 0. 

Remark 2.1 A formal expansion (using a Chapman- Enskog methodology) up to the first 
order in A 2 of the EPB model leads to the following model: 

d t n x + d x (n x u x ) = 0, 
d t {n x u x ) + d x {n x u x u x ) + d x n 

We see that the EPB is a perturbation of the ICE model by a dispersive term with a third 
order derivative in n x . For this reason, in the sequel, the A = 0(1) regime will be referred 
to as the dispersive regime, while the A 1 will be referred to as the hydrodynamic regime. 



-\ 2 n x d x (—d 2 x (logn ; 



0(A 4 ). 
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3 Time semi-discretization, AP property and linearized 
stability 



3.1 Time-semi-discretization and AP property 

We denote by 5 the time step. For any function g(x,t), we denote by g m (x) an ap- 
proximation of g(x,t m ) with t m = mS. We present two time-semi-discretizations of the 
problem. The first one is based on the EPB formulation, and the second one, on the 
REPB formulation. 



3.1.1 Time-semi-discretization based on the EPB formulation 

Classically, when dealing with Euler-Poisson problems, the force term in the momentum 
equation is taken implicitly. In the case of the two-fluid model (when the electrons are 
modeled by the compressible Euler equations instead of being described by the Boltzmann 
relation), S. Fabre has shown that this implicitness is needed for the stability of the scheme 
(an explicit treatment of the force term leads to an unconditionally unstable scheme |18j). 
Additionally, this implicitness still gives rise to an explicit resolution, since the mass 
conservation can be used to update the density, then the Poisson equation is used to 
update the potential, and finally the resulting potential is inserted in the momentum 
equation to update the velocity. 

We will reproduce this strategy here and consider the following time-semi-discretization 
based on the EPB formulation: 

<TV A ' m+1 - n A ' m ) + d x (n x ' m u X ' m ) = 0, 

cT 1 (72 A ' m+ V' m+1 - n x ' m u x ' m ) + d x (n x ' m u x ' m u x ' m ) = n x ' m+1 d x (f) X ' m+1 , 

A 2 9 2^+l = n X,m+l _ e -4><«+\ 

This scheme is Asymptotic-Preserving. Indeed, letting A — > in this scheme with a 
fixed 5 leads to 

< T 1 (n°' m+1 - n°' m ) + d x {n ' m u ' rn ) = 0, 

<r i( n o 1 m+i u o 1 m+i _ n o,m M o,m) + d x (n°' m u 0,m u°' m ) = n°' m+1 d x (j)°> m+1 , 



= n°' n+1 -e-* o ' m+1 . (3.1) 

By using (13.1 j) and the same algebra as for (I2.19p . we find that this scheme is equivalent 
to 

<T 1 (n°' m+1 - n°' m ) + d x {n ' m u°' m ) = 0, 

6 -l^ n 0,m+l u 0,m+l _ n 0,m M 0,m^ + ^O.^O.^O.m) + ^(n°' m+1 ) = 0, 

which provides a semi-implicit discretization of the ICE model, with an implicit treatment 
of the pressure term in the momentum conservation equation. 
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However, when the scheme is discretized in space, the algebra leading to (12.191) is 
no longer exact. Let us denote by _D0 o,m+1 the discretization of the space derivative 
operator d x (jP ,m+1 . Then, the limit A — > of the fully discrete scheme gives rise to the 
approximation n 0,m+1 -D(hi?7 ' m+1 ) of the space derivative d x (n 0,m+1 ) instead of the natural 
derivative Dn°' m+1 . In particular, this expression is not in conservative form. Therefore, 
the use of this scheme may lead to a wrong shock speed if shock waves are present in the 
solution. 

Another drawback of this scheme is the lack of pressure term in the momentum equa- 
tion. As a consequence, the hydrodynamic part of the model is a pressureless gas dynamics 
model, which is a weakly ill-posed model (with, e.g. the possibility of forming delta con- 
centrations [2j [3j IU El [7] ) . The weak instability may lead to spurious oscillations in the 
solution. 

For these reasons, another scheme, based on the REPB formulation, is proposed in 
the next section. 



3.1.2 Time-semi-discretization based on the REPB formulation 

We reproduce the same strategy (i.e. an implicit evaluation of the force term in the 
momentum conservation equation) starting from the REPB formulation. This leads to 
the following scheme: 

5- 1 ((n x ' rn+1 u x ' m+1 ) - (n A '"V' m )) + d x (n x ' rn u x ' m u x ' m ) + d x n x ' rn = 

= \ 2 d x (d 2 x <p x > m+l + l -{d x <t> x > m+1 )^j , (3.3) 

Formally passing to the limit A — )• with fixed 5 in this scheme leads to the following 
scheme: 

( T 1 ((ri > m+ V' m+1 ) - (n°' m M °' m )) + d x (n > m u°> m u°' m ) + d x n°' m = 0, (3.6) 
= n°' m+1 -e-^ m+1 . (3.7) 



Eqs. (13. 5p . (13. 6p are the standard time-semi-discretization of the ICE model. We now 
note that the pressure term d x n°' m is explicit (it was implicit in the scheme based on the 
EPB formulation). Additionally, if a space discretization is used, the discretization of this 
term will stay in conservative form, by contrast to the EPB-based scheme. Finally, the 
hydrodynamic part of the scheme (I3.2p - (I3.4I) is based on the ICE model, not on the pres- 
sureless gas dynamics model. Therefore, its discretization will avoid the possible spurious 
oscillations that might appear in the EPB-based scheme in the presence of discontinuities 
or sharp gradients. 
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3.2 Linearized stability analysis of the time-semi-discretization 

The goal of this section is to analyze the linearized stability properties of both schemes. 
More precisely, we want to show that both schemes are stable under the CFL condition 
of the ICE model, irrespective of the value of A when A — > 0. This property is known 
as 'Asymptotic-Stability' and is a component of the Asymptotic-Preserving property (see 
section [1]). Indeed, the faculty of letting A — > in the scheme with fixed 5 is possible only 
if the stability condition of the scheme is independent of A in this limit. We will prove 
L 2 -stability uniformly with respect to A for the linearization of the problem f l2.14p - fl2.16p . 

In general, time semi-discretizations of hyperbolic problems are unconditionally un- 
stable. This is easily verified on the discretization (13. 5ft - (13. 7p of the ICE model. This is 
because the skew adjoint operator d x has the same effect as a centered space-differencing. 
For fully discrete schemes, stability is obtained at the price of adding numerical viscosity. 
To mimic the effect of this viscosity, in the present section, we will consider the linearized 
Viscous Euler-Poisson-Boltzmann (VEPB) model, which consists of the linearized EPB 
model f l2.1ip -f f2T3|) with additional viscosity terms (in this section, we drop the tildes for 
notational convenience): 

d t n x + d x u x - fid 2 x n x = 0, 
d t u x - (3d 2 x u x = «9 X A , 
A 2 ^ A = n x + A . 

where /3 is a numerical viscosity coefficient. We keep in mind that, in the spatially 
discretized case, (3 proportional to the mesh size h: 

(3 = ch, (3.8) 

with the constant c to be specified later on. Similarly, the linearized Reformulated Viscous 
Euler-Poisson-Boltzmann (RVEPB) model is written: 

d t n x + d x u x - Pd 2 x n x = 0, 
d t u x + d x n x - fid 2 x u x = X 2 d 3 x (f) x , 
\ 2 d 2 x <p x = n x + 4> X - 

The time discretization of these two formulations (which are also linearizations of the 
EPB or REPB-based schemes with added viscosity terms) are given by 



<T 1 (n A ' m+1 - n x ' m ) + 3 x u x ' m - (3d 2 x n x ' m = 0, (3.9) 

A 2 <9 2 A ' m+1 = n x ' m+1 + A ' m+1 , (3.11) 
for the EPB-based scheme and by 

r i(y, m +i - n x ' m ) + d x u x ' rn - (3d 2 x n x ' m = 0, (3.12) 

S -l( u \,m+l _ M A,m) + dxU X,m _ ^^A.m = X 2 d^ X ' m+1 , (3.13) 

X 2 d 2 x (j) x ' m+1 = n x ' m+1 + A - m+1 , (3.14) 
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for the REPB-based one. 

Passing to Fourier space with £ being the dual variable to x, and eliminating A,m+1 , 
we find the following recursion relations: 

for the EPB-based scheme and 

1 + A 2 4 

for the REPB-based one. 

The characteristic equations for these two recursion formulas are 



q 2 - 2q(l - Pet ~ n J,\m ) + (1 - /W = 0, (3.15) 



^2 

2(1 + A 2 f 2 
and 



q 2 - 2q(l - P^ 2 S + 2(1 A ^ 2) ) + (1 - + = 0, (3-16) 

respectively, where q is the characteristic root. Each of these quadratic equations has two 
roots q±(£) which provide the two independent solutions of the corresponding recursion 
formulas. Their most general solution is of the form 

(S3l!)-?<*or(» 

where n±(£) and depend on the initial condition only. A necessary and sufficient 

condition for L 2 stability is that |<?±(£)| < 1. However, requesting this condition for all 
£ G R is too restrictive. To account for the effect of a spatial discretization in this analysis, 
we must restrict the range of admissible Fourier wave- vectors f to the interval [— ? , £]. 
Indeed, a space discretization of step /t cannot represent wave- vectors of magnitude larger 
than j-. This motivates the following definition of stability: 

Definition 1 The scheme is stable if and only if 

|g±(OI < h Vf such that |f | < |. (3.17) 



Now, our goal is to find sufficient conditions on 5 such that either schemes are stable. 
More precisely, we prove: 



14 



Proposition 2 For both the EPB-based scheme l\3.9\) - /[3~Tl\) or the REPB-based scheme 
A3.1fy) - [3.14\ )> there exists a constant C > independent of X when A — > such that if 



5 < Ch, the scheme is stable. 

This condition states that the schemes are stable irrespective of how small A is. We 
say that the schemes are 'Asymptotically-Stable' in the limit A — > 0. We note that this 
stability condition is similar to the CFL condition of the ICE model, which is the limit 
model when A — > 0. 

Proof: We first define conditions such that the constant term of the quadratic equations 
( 13.151) or ( I3.16P is less than 1. For the EPB-based scheme (see 13.151) . this condition is 
5 < 1 /Pi 2 , for all f such that |f| < n/h. For reasons which will become clear below, we 
rather impose: 



7T 



S<^, Ve such that |e|<-, (3.18) 



which, with (13.81) . is equivalent to: 



5 < C x h, with d = —r. (3.19) 

ZC7T 2 



For the REPB-based scheme (see 13.161) . this condition is 



or, with CT) : 



5<^|^, V£ such that |e|<^, 



5 < Cxh, with d = — (3.20) 

1 + C 2 7T 2 



Now, under these conditions, we are guaranteed that the two roots satisfy ( I3.17P if and 
only if the discriminant of the quadratic equation is non-positive. Indeed, in this condition, 
the two roots are conjugate complex numbers and their product, i.e. the square of their 
module, which is equal to the constant term of the quadratic equation, is less than one. 

It is a matter of computation to check that the discriminant has the same sign as the 
expression F(5) given by: 

F(5) =5 2 + 4/3(1 + X 2 e)S - 4 1+ ^ 2 , (3.21) 

in the case of the EPB-based scheme and by 

F(6) = 6 -4(3 A2 ^ 2 6-4 ^ , 

in the case of the REPB-based scheme. 
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In the case of the EPB-based scheme, we use (13 . 1 8[) to estimate the second term of 
F(5) in dS2H): 



1 -i- \ 2 f= 2 
F(5)<5 2 - 2 i±^, 

and a sufficient condition for F(5) to be non- positive is that § < y/2(l + A 2 ^) 1 / 2 ^ -1 . 
This relation is true for all £ such that |£| < Tr/h if 6 < v / 27r~ 1 (/?, 2 + A 2 ^ 2 ) 1 / 2 and a 
sufficient condition is that 5 < C 2 h, with C 2 = \f2ix~ l . Now, taking C = min{Ci,C2} 
with C\ given by (13 . 19[) leads to the result. In fact, the optimal numerical viscosity is 
such that d = C 2 , i.e. c = (2y/2n)- 1 . In the case of the REPB scheme, we estimate 
F{5) by 



1 4- \ 2 F 2 



and a sufficient condition for F(5) to be non- positive is that 5 < 4/3(1 + A 2 £ 2 )(A 2 £ 2 ) -1 . 
In view of (13.81) . this relation is true for all £ such that |£| < n/h if 5 < Ac(\ 2 ir 2 )~ l h{h 2 + 
\ 2 7[ 2 yi 2 and a sufficient condition is that 5 < C 2 h, with C 2 = 4c. Now, taking C = 
min{Ci,C2} = C\ with C\ given by (I3.20p leads to the result. The optimal numerical 
viscosity can be chosen to minimize C, which leads to c = 7r _1 . This ends the proof of 
statement El ■ 

As a conclusion, we can see that both the EPB-based and REPB-based schemes have 
similar Asymptotic-Stability properties as A — > 0. Therefore, they must be selected 
on the basis of other criteria. The fact that the REPB-based scheme has a well-posed 
hydrodynamic part and leads to a discretization of the ICE model in conservative form 
in the limit A — > are indications that this scheme should be preferred to the EPB-based 
scheme. In the next section, we will present numerical results that support this statement. 



4 Spatial discretization 

We introduce (Cj)^ =1 a uniform subdivision of the computational domain Q G R such 
that Q = UjLiCj. The interface between Cj and C J+ i is the point Xj+i/2- We denote by 
U™ the approximate vector of the density and momentum at time t m on the cell Cj, 




We use a time-splitting method to compute the density and momentum at time t 
There are three steps to pass from U m to JJ m+1 which are described below. 
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4.1 Hydrodynamic part 

The first step of the splitting is the finite- volume computation of the state such that 

rr# _ Tjm pm _ pm 

U U j , r 3+1/2 J j-l/2 _ n 

5 + h 

where i^+1/2 is the numerical flux computed at the interface Xj+1/2. 

We have used a Local Lax-Friedrichs [22] (or Rusanov [28] or degree polynomial 
[16] ) solver. This solver is an improved version of the Lax-Friedrich solver which has been 
successfully used in conjunction with AP-schemes for the two-fluid Euler-Poisson problem 
(see [11] ) or for small Mach- number flows [17]. This solver depends on a local estimate of 
the maximal characteristic speed. This estimate proceeds as follows. We introduce 

(a% 1/2 = max (u™ 1/2 + 1, uf +1 + l) and (a~)f +1/2 = min (u™ - 1, uJ +1/2 - l) , 

where u™ +1 / 2 = (u 1 J l + uJ L +1 )/2. Then, the local maximal characteristic velocity is estimated 
by 

aJ +1/2 = max (\(a~)™ 1/2 |, |(a+)£ 1/2 |) , 
and the numerical flux at 2^+1/2 is given by: 

F j+ i/2 = \ {F{Uf) + F{Uf +1 ) + af +l/2 (U™ - Uf +1 )) . (4.1) 

The time step 5 must satisfy the CFL condition | < max a™ +1 i 2 to ensure stability. In 
practice, the time step is chosen at each iteration to enforce this stability condition. As 
for boundary conditions, we impose fictitious states Ui and U r across the left and right 
boundaries respectively and compute the corresponding fluxes across the boundaries using 
the same formula (14. ip . 

4.2 Potential update 

The second step of the time-splitting is the computation of the potential. We use n* to 
compute <p m+1 with the discretized Poisson equation given by a finite difference approxi- 
mation of the Poisson equation. 

_ + ™+i) + e-^ +1 = n f. 

The boundary conditions are given by <\>i = — log nf on the left hand side of the domain 
and 4> r = — log nf on the right hand side of the domain. The non- linear system is solved 
with newton method. This iterative algorithm needs a good initial guess of the solution 
to be efficient. The initial guess is the potential at previous time <j) m . For the first step, 
we choose the quasi-neutral potential (0°) = (— logn°) as an initial guess. 
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4.3 Source term 



In the REPB form, the source term Q at the right-hand side of H3 .31) can be written 
according to: 

Q = X 2 {<)A<to) I »,o ;)?:>)■ 

This expression is discretized thanks to a finite difference approximation. For a cell Cj in 
the domain, the source term Qj is given by centered finite difference approximation: 

A 2 

Qj = 7^3 (Oj+2 - 20 i+ i + 20j_i - (pj- 2 ) + (<f>j+i ~ 2<Pj + 4>j-i) (<f>j+i ~ <t>j-i)) ■ 

On the first cell C\ the source term is computed using a decentered finite difference 
approximation: 

A 2 / 1 

Qi = ^ ( (03 - 30 2 + 30! -00+2 - 2 ^ + ^) - ^ 

and similarly in the last cell. 



5 Numerical results 

We present three classes of numerical results : the first test case is a solitary wave travelling 
in a plasma. This test case shows the ability of the numerical schemes to handle the 
dispersive regime. The second test case is related to the quasi-neutral limit A — > of 
the Euler-Poisson-Boltzmann model: it is a Riemann problem to check the ability of the 
scheme to handle hydrodynamic phenomena like shocks. The last test case has been 
previously investigated by Liu and Wang in [231 121] and corresponds to the occurrence of 
multi-valued solutions in the semi-classical setting. 

5.1 Soliton test case 
5.1.1 Description 

The dispersive nature of the Euler-Poisson system is shown in [25]. Therefore, the EPB 
system, like other nonlinear dispersive models such as the KdV equation, exhibits solitary 
wave solutions. These special solutions are particularly convenient to test the ability of 
the EPB and REPB schemes to capture the dispersive regime. Solitary waves also provide 
interesting quantitative checks. Indeed, while travelling through the plasma, the soliton 
maintains its shape and velocity. Therefore, one can check the accuracy of the numerical 
schemes by observing how well they preserve the soliton shape and velocity over time. 

We now summarize the establishment of this special solution. We refer to [5] for a de- 
tailed description and physical considerations. For this derivation, we use non-dimensional 
units and we now precise the corresponding scaling units. The space scale related to these 
solitons is the Debye length Ad. For this reason, we take A = 1 in all this section. The 
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size of the computational domain is equal to several Debye lengths (about 50Xd are used 
in the subsequent simulations). The density of the undisturbed plasma (away from the 
support of the solitary wave) is chosen as the characteristic density and in dimensionless 
units, is equal to 1, so that the electron density is equal to where is the electrostatic 
potential energy. 

In the frame moving with the wave, we denote by n s ,u s , <p s the density, velocity and 
potential of the plasma. These quantities are constant in time, and satisfy the following 
relations: 

d x (n s u s ) = 0, 
d x (n s u s u s ) = n s d x (j) s , 
d 2 x <p s = n s -e-^. 

The momentum being uniform in x, we write q = n s u s = UoUq where uq = n s (0) and 
no = u s (0). The momentum conservation law can be written as d x q 2 /n s = n s d x cf) s . 
Consequently, we have d x (q 2 /n 2 ) = 2d x (p s . For all x G [0,x max ], we get: 



2 \n 2 J v ; 2 \n 2 s/ 

The potential being defined up to an additive constant, we choose this constant such that 
0s(O) = 0. The ion density is then given by 

-.M=Gr!ir /2 > (5 - 2) 

In the present analysis, we assume that no = 1, i-e. no is equal to the density of the 
undisturbed plasma. Inserting this relation in the Poisson equation yields the following 
equation for the potential: 

d 2 x <Ps=(l + 2 4) V2 -e-*>, (5.3) 

with the condition 

0,(0) = 0. (5.4) 

One needs another condition at x = to set up a Cauchy problem for (15.31) . Note that 
<p s = is an obvious solution of equation (15. 3p and satisfies the homogeneous condition 
<9x0s(O) = 0. It corresponds to the state of the undisturbed plasma. To capture a non- 
trivial solution, we must alter this condition by a small disturbance, setting it to 

d x( j) s (0) = 77, (5.5) 

with 77 'small'. 

The behavior of these solutions can be clarified thanks to the theory of the Sagdeev 
potential (which is a primitive with respect to <ft s of the right-hand side of (15. 3p ). Details 



19 



on this study can be found in [5]. Here we just recall that shock waves in a cold-ion 
plasma can exist only for 1 < uq < 1.6 (Bohm criterion). The sign of r\ determines the 
type of solution which can be found: a positive rj leads to potential barrier that forms a 
sheath, whereas a negative value gives rise to a monotonic transition to a negative <fi which 
forms a solitary wave corresponding to a potential and density disturbance propagating to 
the right (for instance) with velocity uq. There is no analytic solution to equation (15. 3p . 
Resorting to numerical simulation is the only way to determine these solutions. Details 
about this numerical method are given below. Once the soliton potential is known, the 
density is computed thanks to (15.21) . Since the present analysis has been performed in a 
co-moving frame with the soliton, moving with velocity uq, the velocity in the laboratory 
frame is u s + uq with u s = q/n s . In the subsequent numerical simulations n s , uo + u s , <fi s 
are taken as an initial condition for the scheme. 

5.1.2 Numerical results for the soliton test case 

In this test case, where the Debye length and the size of the computational domain are of 
the same order, both the EPB and REPB schemes are expected to be correct. However, 
this test case provides a way to achieve quantitative comparisons of the numerical solutions 
to an analytical reference solution. These comparisons permit to quantify the order of 
accuracy and amount of numerical diffusion of the two schemes. The analytical solution 
is easily obtained at time t by a spatial translation of the solution at time of a distance 
u t. 

This test case is implemented as follows. The boundary conditions are periodic, which 
ensures that the soliton can be tracked on long simulation times without the need for a 
large computational domain. The length L of the computational domain is defined in 
relation to the choice of the initial condition as explained below. We denote by = L/uq 
the travel time of the soliton in the domain [0, L] . First a numerical solution of (15.31) 
with initial conditions (15 .4p . (15. 5p is computed with an explicit finite difference scheme on 
a mesh of step Ax rei and provides the reference solution. Ax rei is chosen small enough to 
provide an almost 'exact' reference solution. 

For suitable r\ and u given by the study of the Sagdeev potential [5] , the potential 
oscillates in space for positive x. One wants a single potential well to initialize the scheme. 
To this aim the number of nodes N is taken such that the initial condition shows a single 
peak. Moreover, it is such that <f) T jf iel is close enough to to ensure that periodic boundary 
conditions will be accurate enough. This reference solution is interpolated to provide an 
initial condition for the EPB and REPB schemes, and to compute numerical errors on 
the density, momentum and potential. 

The results of this comparisons are now commented. Figure [1] and [2] show the density 
and velocity in a soliton at times 0, £l/5, 2^/5 and ti, computed with the reformu- 
lated scheme. The shape of the initial condition is conserved with time, but the peak is 
damped and does not return to its original location after t L . This effect is due to the 
numerical diffusion, which is inherent to the numerical method, and can be observed with 
the classical scheme as well. In order to accurately compare the reformulated and classical 
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schemes, one needs to perform a more precise study of this damping. The forthcoming 
convergence study uses six grids, with a range of space steps from h to h/64, where h 
corresponds to 250 cells. It confirms that both schemes are first order in space, and even 
if the REPB scheme suffers from a larger numerical diffusion than the classical one, both 
provide satisfactory results for this test case. 



Density in a soliton Momentum in a soliton 




°' 8 L/5 2L/5 3L75 4L/5 L °' 2 L/5 2L75 3L/5 4L/5 L 

Position Position 

Figure 1: Density (left) and momentum (right) as a function of space, for a soliton moving 
to the right, at time 0, £l/5, 2t^/5 and t^, computed with the REPB scheme. At time 
the wave has travelled through the whole computational domain. The damping of the 
density and velocity peaks are clearly visible due to the large space step, Ax = h/4. 

The numerical convergence of the schemes in space is investigated by comparison with 
the reference solution; at time t/5 the L°° relative error with the reference solution is 
computed. For instance the density error is 



max 
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where n re f is the density of the reference solution, i.e. the initial density translated by L/5, 
and n num is the density computed with the classical scheme. Such errors are defined for 
the reformulated scheme and for the momentum and potential. The schemes are tested 
on six grids, made of 250, 500, 1000, 2000, 4000, 8000 and 16000 cells. Table [Hand figure [3] 
confirm that both numerical schemes are first order in space. The REPB scheme suffers 
from a larger numerical dissipation than the EPB scheme. 

The numerical dissipation of the schemes can be measured in another way. Indeed, the 
amplitude of the numerical soliton is damped with time. The following study quantifies 
the damping rates of both the EPB and REPB schemes. This study is performed on 
a long time simulation: its final time is 2t^. This study compares the damping rate of 
the density amplitude over one time increment At (not related with the numerical time 
step), at different times of the simulation. The time increment At is £l/5. The density 
amplitude n max ((/c + l)At) is compared to the density amplitude at the previous time 
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Density in a soliton Momentum in a soliton 
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Figure 2: Density (left) and momentum (right) as a function of space, for a soliton moving 
to the right, at time 0, £l/5, 2£l/5 and t^, computed with the REPB scheme. At time tz, 
the wave has travelled through the whole computational domain. The solution at time t E 
can be superimposed to the solution at time 0, since the grid is fine enough (Ax = h/64) 
and the numerical diffusion is very small. 



N 


£EPB(n) 


£REPB(n) 


e E PB(nu) 


£REPB(nu) 


£epb(0) 


£repb(4>) 


h 


0.053 


0.102 




0.111 


0.215 


0.066 


0.116 


h/2 


0.028 


0.060 




0.060 


0.131 


0.028 


0.066 


h/A 


0.014 


0.035 




0.032 


0.073 


0.014 


0.035 


h/8 


7.6 x 10~ 3 


18.5 x 10' 


-3 


1.64 x 10^ 2 


3.87 x 10~ 2 


7.12 x 10~ 3 


18.4 x 10~ 3 


h/lQ 


3.86 x 10~ 3 


9.57 x 10' 


-3 


8.40 x 10~ 3 


2.00 x 10~ 2 


3.50 x 10~ 3 


9.42 x 10~ 3 


h/32 


1.76 x 10" 3 


4.67 x 10' 


-3 


3.83 x 10~ 3 


1.00 x 10" 2 


2.31 x 10~ 3 


4.71 x 10" 3 


h/QA 


8.89 x 10" 4 


2.37 x 10" 


-3 


1.97 x 10" 3 


5.09 x 10" 3 


1.07 x 10" 3 


2.41 x 10" 3 



Table 1: Comparison of the error in L°° norm for the density, momentum and potential 
with various space grid sizes using the EPB and REPB schemes. The computation of the 
error is made with the analytical solution for the soliton test case 
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Relative error on density 



Relative error on momentum 




Ax 

Relative error on potential 




Figure 3: Comparison of the error in L°° 
norm for the density, momentum and po- 
tential with various spacial mesh sizes Ax 
for the EPB and REPB schemes. Both 
schemes are first order in space. 
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increment n max (&At) for < k < 9. We measure the decrement of the amplitude A& 
between kAt and (k + l)At as follows: 

(5.6) 

The damping rate of the wave amplitude appears on figure H] for both the EPB and REPB 
schemes. Two spatial grids are used for this comparison: a coarse grid with 1000 cells 
and a fine grid with 16000 cells. The REPB scheme shows a larger dissipation than the 
classical one. The evolution of the damping rates with simulation time is similar for 
the two schemes : on a coarse grid the peak is damped faster at the beginning of the 
simulation. 

Damping on soliton (density amplitude ) 
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Figure 4: Comparison of the damping rates for the density amplitude of the soliton for 
the EPB and REPB schemes on two different spatial mesh sizes: a coarser grid with 1000 
cells and a finer grid with 16000 cells are used for both schemes. 
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5.2 Riemann problem 

The previous test case, in which the domain size and Debye length are of the same order 
of magnitude, was designed for the dispersive regime A = 0(1). The present test case 
aims at investigating how the schemes perform in the hydrodynamic (or quasi-neutral) 
regime A 1. Therefore, in this test case, values of A ranging from small to very small 
small are used. When A is small, the REPB model explicitly appears as a perturbation 
of the ICE model, which is a strongly hyperbolic and conservative model. Therefore the 
REPB-based scheme should be accurate in the hydrodynamic regime, in particular for the 
computation of solutions involving discontinuities. By contrast, the hydrodynamic part 
of the EPB model is a pressureless gas dynamics model, which is not strictly hyperbolic, 
and is thus weakly unstable. For this reason, the EPB-based scheme is expected to be 
less accurate in the small Debye length limit. The present test problem aims at testing 
the validity of these predictions. 
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The test case is a shock tube problem involving two outgoing shock waves. The initial 
density is a constant, while the initial velocity has a jump at x = between the constants 
ul = +1 and ur = — 1. The density in the intermediate state of the Riemann problem 
is larger and the velocity is zero. Two outgoing shock waves appear on each side of this 
intermediate state. The computational domain is [—0.2; 0.2]. The dimensionless Debye 
length A varies from 10~ 2 to 10~ 4 . 

The value A = 10 -2 is large enough for singularities near the shock waves to appear, 
due to the dispersive nature of the Euler-Poisson-Boltzmann model. In the semi-classical 
setting, the framework of multi- valued solutions [231 121] can be used to explain the qual- 
itative features of the classical solutions. Indeed, classical solutions keep a signature of 
these underlying multi-valued solutions, in the form of singularities (when the number 
of branches changes) and oscillations (when several branches co-exist and the solution 
'oscillates' between these branches). In the A — > limit, the entropic solutions of the 
ICE model capture the average value of these oscillations, but not the details of them. 
As we will see, the EPB-based scheme keeps better track of these oscillations but when 
the mesh size does not resolve the Debye length, the oscillations become mesh-dependent. 
On the other hand, the REPB-based scheme directly provides the entropic solution of the 
limiting ICE model and is better suited to capture the average value of these oscillations 
(i.e. the weak limit of the solutions of the EPB model when A — > 0). 

Fig. [5] confirms that the numerical solution provided by the EPB-based scheme in 
under- resolved situation is mesh-dependent. It displays the density computed by the 
EPB-based schemes for two different mesh sizes. The density peaks computed on the finer 
mesh are much finer and higher than those computed on the coarse mesh. By contrast, 
the velocity and potential remain finite regardless of the mesh size (not displayed). Fig. 
[6] shows the velocity and potential computed with the EPB-based scheme compared to 
the REPB-based scheme. Like in the soliton test case, the REPB-based scheme shows a 
slightly larger numerical dissipation. 

The value A = 10 -4 is small enough to observe the hydrodynamic regime. Both 
schemes show an accurate determination of the shock speed but the EPB scheme leads 
to spurious oscillations in the neighborhood of the shock. These oscillations are shown 
on figure [7J They are mesh dependent and occur even with smaller dimensionless Debye 
length A. Whatever the mesh size is the solution computed with the reformulated REPB 
scheme does not present such oscillations. 

This test case shows that the EPB and REPB-based schemes have a very different 
behavior when A 1. In such under-resolved situations, while the EPB scheme presents 
mesh-dependent oscillations of finite amplitude, the REPB-based scheme provides an ac- 
curate approximation of the entropic solution of the limiting ICE model. As a conclusion, 
the REPB-based scheme should be preferred for the A < 1 regime. 
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Riemann Shock Waves, t = 0.2, X = 10 



Riemann Shock Waves, t = 0.2, % = 10 
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Figure 5: Density as a function of space for the Riemann shock-wave problem at time 
t = 0.2 with dimensionless Debye length A = 10~ 2 . Left: classical scheme with a coarse 
grid (N = 2000). Right: classical scheme with a fine grid (N = 32000). Density peaks 
are sharpened when the mesh size decreases 
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Figure 6: Velocity (left) and potential (right) as a function of space for the Riemann shock- 
wave problem at time t = 0.2 with dimensionless Debye length A = 10 -2 . EPB-based 
(dashed line) and REPB-based (solid line) schemes on a fine grid (N = 32000). 
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Riemann Shock Waves, time t=0.1 and \=0.2,X = 10"' 



Riemann Shock Waves, time t=0.1 and t=0.2, ^ = 10' 
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Figure 7: Density as a function of space for the Riemann shock-wave problem at time t = 
0.1 and t = 0.2 with dimensionless Debye length A = 10~ 4 . Left: EPB and REPB-based 
schemes on a coarse grid (N = 2000). Right: same schemes on a fine grid (N = 32000). 



5.3 Dispersive solutions test-cases 
5.3.1 Description 

In this section, we present test-cases which are inspired from [24]. In [24], the goal was to 
explore the computation of multi-valued solutions in the semi-classical setting by means 
of the level-set method. Here we consider only classical solutions. The first test case is 
referred to as a five branch test case (because it corresponds to the occurrence of a five 
branch multi- valued solutions in the semi-classical setting). The second test is a seven 
branch test case (for the same reason). The initial densities of both test-cases are bumps 
with a gaussian shape. This shape leads to a potential well which generates an induced 
electric field. This electric field in turn contracts the density bump leading to a positive 
feedback amplification. This effect can be further amplified by setting up appropriate 
initial velocities. In the dispersive regime the amplification of the density peak can lead 
to singularities, whereas in the hydrodynamic regime the density spreads out in the entire 
computational domain and its profile remains smooth. 

The emergence of singularities has been investigated by Liu and Wang [2"3"| |2"4] . The 
authors compute multi- valued solution for similar test cases but with the standard Poisson 
equation (without the exponential term coming from the Boltzmann relation) that allows 
the computation of analytical solutions. Here, no analytical solution is available but the 
two schemes (EPB and REPB) are tested one against each other and against a numerically 
computed reference solution on a very fine mesh. Due to the singularities appearing when 
A = 1 the numerical error is computed with the potential, which remains finite in every 
situation. The tests are run with A = 1 and A = 10~ 2 in order to explore both the 
dispersive and hydrodynamic regimes. The results would be similar if A was further 
reduced. 
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5.3.2 Five-branch solution 



The initial condition for this first test-case is given by 



n 
u 



7T 

sin 3 x. 



The computational domain is [0, 27r]. The initial density and velocity appear on figure 
[8] and [91 together with the numerical solutions computed at time t = 1 by means of the 
EPB and REPB-based schemes. Table [2] shows the numerical error on the potential by 
comparison against a solution computed with the classical scheme on a fine grid. 



5 branches test case, time t = 1 , X = 10 



5 branches test case, time t = 1 , X = 10 
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Figure 8: 
test case, 
computed 



Density (left) and velocity (right) as a function of space for the five-branch 
The dimensionless Debye length A is 1. Numerical solutions at time t — 1 are 
with the EPB and REPB-based schemes on a grid with 2000 cells. 



5 branches test case, time t = 1 , X = 1 ' 



5 branches test case, time t = 1, X = 10 ' 
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Figure 9: Density (left) and velocity (right) as a function of space for the five-branch test 
case. The dimensionless Debye length A is 10 -2 . Numerical solutions at time t = 1 are 
computed with the EPB and REPB-based schemes on a grid with 2000 cells. 
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N 


£epb(<P), A = 1 


£repb((^), A = 1 


Sepb(<P)A = 10~ 2 


£repb{4>), A = 10~ 2 


8000 


5 x 10" 5 


1.4 x 10" 4 


4.0 x 10" 4 


7.6 x 10" 4 


4000 


1.5 x 10" 4 


2.8 x 10" 4 


1.0 x 10" 3 


1.8 x 10~ 3 


2000 


3.4 x 10~ 4 


5.4 x 10~ 3 


2.3 x 10~ 3 


3.7 x 10~ 3 



Table 2: Five-branch test case : comparison of the error in L°° norm on the potential 
for various grid sizes and dimensionless Debye length A using the EPB and REPB-based 
schemes. 

5.3.3 Seven-branch solution 

The initial condition for this second test-case is given by 

u = sin(2x) cos x. 

The initial density and velocity appear on figure dD] and [TU together with the numerical 
solutions computed at time t — 1 using the EPB and REPB-based schemes. Table [3] 
shows the numerical error on the potential by comparison against a solution computed 
with the classical scheme on a fine grid. 



N 


£epb(0), A = 1 


£repb{4>), A = 1 


Sepb{4>)A = KT 2 


Srepb{4>), A = 10~ 2 


8000 


4.6 x 10" 5 


4.5 x 10" 4 


2.3 x 10~ 4 


7.6 x 10" 4 


4000 


1.4 x 10~ 4 


7.4 x 10" 4 


6.7 x 10~ 4 


1.6 x 10~ 3 


2000 


3.2 x 10~ 4 


1.2 x 10~ 3 


1.29 x 10" 3 


3.5 x 10~ 3 



Table 3: Seven-branch test case : comparison of the error in L°° norm on the potential for 
various grid sizes and dimensionless Debye length A using the classical and reformulated 
Euler-Poisson-Boltzmann schemes. 



5.3.4 Analysis of the results for the five and seven branch test cases 

In both the dispersive (A = 1) and hydrodynamic (A = 10~ 2 ) regimes, the two schemes 
give similar results. On figure [SI [HI QUI and [TU overlapping lines for the solution at time 
t — 1 computed with the EPB and REPB-based schemes confirm this similar behavior. 
One can exhibit some differences thanks to the numerical convergence study. Tables [2] 
and [3] show the same differences between the two schemes as in the previously discussed 
soliton test case. The error of the reformulated scheme is slightly larger than that of the 
classical scheme, and this difference is more obvious when A = 1. 

6 Conclusion 

In this paper, we have analyzed two schemes for the Euler-Poisson-Boltzmann (EPB) 
model of plasma physics, and compared them in different regimes characterized by dif- 
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7 branches test case, time t = 1 , X - 1 7 branches test case, time t = 1 , X - 1 0" 
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Figure 10: Density (left) and velocity (right) as a function of space for the seven-branch 
test case. The dimensionless Debye length A is 1. Numerical solutions at time t — 1 are 
computed with EPB and REPB-based schemes on a grid with 2000 cells. 
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Figure 11: Density (left) and velocity (right) as a function of space for the seven-branch 
test case. The dimensionless Debye length A is 10~ 2 . Numerical solutions at time t = 1 
are computed with EPB and REPB-based schemes on a grid with 2000 cells. 
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ferent values of the dimensionless Debye length A. The dispersive regime corresponds to 
A = 0(1) while the hydrodynamic regime is characterized by A <C 1. When A — > 0, the 
EPB model formally converges to the Isothermal Compressible Euler (ICE) model. The 
first scheme we have considered is based on the original EPB formulation of the model. 
The second one uses a reformulation (referred to as the REPB model) in which the model 
more explicitly appears as a singular perturbation of the ICE Model. 

We have provided a stability analysis of the two schemes, showing that both schemes 
are stable in both the dispersive and hydrodynamic regimes, with stability constraints 
on the time and mesh steps which are independent of A when A — > 0. Finally, we have 
tested them on three different one- dimensional test problems. The first test problem, the 
soliton test, provides an analytical solution in the dispersive regime. The second test 
problem, the Riemann problem with two expanding shock waves, is suitable to explore 
the hydrodynamic regime. Finally, the third test problem involves singularity formation 
in the dispersive regime. 

We have concluded that both scheme have similar behavior in the dispersive regime 
(with a slightly increased, but perfectly acceptable numerical diffusion in the case of the 
REPB-based schemes). By contrast, in the hydrodynamic regime, the EPB-based schemes 
develop oscillations and singularities, which, in under-resolved situations (i.e. when the 
time and space steps are too large to resolve the spatio-temporal variations of the solution) 
prevent any grid convergence of the solution. By contrast, the REPB-based scheme well 
captures the entropic solution of the ICE model, which provides a good approximation of 
the weak limit of the dispersive (oscillatory) solutions of the EPB model in the small A 
regime. 

Future works concern the extension of this analysis to the two- or multi-dimensional 
case, the passage to second order schemes and the pursuit of the analytical investigations 
of the accuracy and stability of the schemes in both regimes. 
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